Fractals and dynamical chaos 
in a random 2D Lorentz gas with sinks 



I. Claus, 1 P. Gaspard, 1 and H. van Beijeren 2 

1 Center for Nonlinear Phenomena and Complex Systems, Faculte des Sciences, 
Universite Libre de Bruxelles, Campus Plaine, Code Postal 231, 
B-1050 Brussels, Belgium 
2 Institute for Theoretical Physics, Utrecht University, 
Leuvenlaan 4, 3584 CE Utrecht, The Netherlands 

Abstract 

Two-dimensional random Lorentz gases with absorbing traps are considered in 
which a moving point particle undergoes elastic collisions on hard disks and an- 
nihilates when reaching a trap. In systems of finite spatial extension, the asymp- 
totic decay of the survival probability is exponential and characterized by an escape 
rate 7, which can be related to the average positive Lyapunov exponent and to the 
dimension of the fractal repeller of the system. For infinite systems, the survival 
probability obeys a stretched exponential law of the form P(c,t) ~ exp(— Ct 1 / 2 ). 
The transition between the two regimes is studied and we show that, for a given 
trap density, the non-integer dimension of the fractal repeller increases with the sys- 
tem size to finally reach the integer dimension of the phase space. Nevertheless, the 
repeller remains fractal. We determine the special scaling properties of this fractal. 
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1 Introduction 



Many recent works have been devoted to the problem of diffusion in the pres- 
ence of static absorbing traps, or sinks. In such models, the diffusing particles 
can be seen as modeling some active chemical species or some kind of excita- 
tion. The trapping then consists in a reaction with another chemical species 
leading to a non-reactive product, or in an excitation trapping converting the 
energy to some other form [1] . The important quantity to study is the survival 
probability P(c, t), i.e., the probability for a diffusing particle to survive in the 
system at time t for a trap concentration c. In the case of randomly distributed 
traps, in a system of infinite spatial extension, the survival probability P(c, t) 
of the diffusing particles obeys asymptotically a stretched exponential law of 
the form [2,3] 



where d is the dimension of the system, c is the concentration of static traps, 
u — — ln(l — c) and is a constant depending on the dimension and on the 
geometry of the system. This behavior is due to the existence with a non- 
vanishing probability of very large trap-free regions in which particles can 
survive for a very long time [3]. 

The stretched exponential law (1) has been much studied in stochastic sys- 
tems although, until now, no study of this phenomenon has been carried out 
in deterministic chaotic systems. Such a study is of special interest in deter- 
ministic chaotic systems where the trapping process is characterized not only 
by the behavior of the survival probability but also by the fractal set of all 
the trajectories that never meet a trap. This set is defined in the phase space 
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of the deterministic system and is referred to as a fractal repeller [4,5]. 

The fractal repeller associated with the deterministic chaotic diffusion in the 
presence of periodically distributed traps or sinks has been studied in a pre- 
vious paper [6]. The system studied there was a two-dimensional periodic 
Lorentz gas: a point particle undergoes elastic collisions on hard disks fixed in 
the plane and forming a regular triangular lattice. Some of these disks, form- 
ing a regular superlattice over the disk lattice, play the role of traps: when the 
point particle collides on one of them, it is annihilated. In this case, since the 
size of trap-free regions is bounded, the decay is exponential, characterized by 
an escape rate 7, so 

P(c,f)~e-*. (2) 

The dynamics is chaotic and the repeller has a non-integer fractal dimension. 

The purpose of the present paper is to introduce a model for the stretched 
exponential law (1) which is based on a deterministic dynamics in order to be 
able to investigate the deterministic features of the process and, especially, the 
properties of the associated fractal repeller. We might suspect that the fractal 
dimension of the repeller is integer and equal to the phase space dimension 
Dh = 3 since, in a random system of infinite spatial extension, the survival 
probability (1) decays more slowly than the exponential (2). However, if the 
dimension is integer, a finer characterization of the repeller is required and we 
might even wonder if the repeller is still a fractal. 

In order to find a solution to this puzzle, we first apply the escape-rate for- 
malism to finite two-dimensional random Lorentz gases. We confirm that, in 
such systems, the asymptotic decay is exponential and that the correspond- 
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ing escape rate can be related to the positive Lyapunov exponent and to the 
fractal dimension of the fractal repeller. 

We next consider infinite 2D random Lorentz gases. Thanks to our fast event- 
driven numerical algorithm and appropriate simulation method, we are able 
to extend the study of the survival probability to these deterministic chaotic 
systems. 

Finally, we study the transition from finite to infinite 2D random Lorentz 
gases, which then allows us to determine the evolution of the scaling prop- 
erties of the fractal repeller when the system size increases and to solve the 
aformentioned puzzle. 

The paper is organized as follows. Section 2 is devoted to finite two-dimensional 
random Lorentz gases with sinks. The asymptotic decay is shown to be ex- 
ponential, described by an escape rate 7. Moreover, we show that this open 
chaotic system is characterized by the existence of a fractal repeller. In these 
open systems, the average Lyapunov exponent and the fractal dimensions of 
the repeller are calculated with a non-equilibrium probability measure. The 
validity of the relation (2) is confirmed for these finite size systems. The case 
of infinite two-dimensional random Lorentz gases is presented in section 3. 
Thanks to our appropriate simulation method, the stretched exponential decay 
is observed for the first time in a deterministic chaotic system. The transition 
from finite systems with exponential decay to infinite ones obeying stretched 
exponential laws is studied in section 4. For a given concentration of traps, the 
Hausdorff dimension of the repeller is shown to increase with system size, to 
finally reach an integer value equal to the dimension of phase space, as it would 
be the case for a closed system. We show that, nevertheless, the repeller re- 
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mains fractal in the infinite-size limit and that it is no longer characterized by 
a simple power-law scaling as the classical fractals but by a more complicated 
scaling. Conclusions are drawn in section 5. 

2 The two-dimensional random Lorentz gas of finite spatial exten- 
sion 

2. 1 Description of the model 

The 2D random Lorentz gas consists of a point particle colliding elastically 
on randomly and uniformly distributed non-overlapping hard disks fixed in 
the plane. In finite (or periodic) geometry we will require in addition that 
the horizon is finite, i.e., that the point particle, starting from any initial 
condition, collides on a disk after a finite time. Under these conditions the 
system is diffusive, with a well-defined diffusion coefficient. Some disks are 
chosen to be sinks: when the point particle collides on one of them, it is 
absorbed, according to the following reaction scheme: 

X + inert disk <-> X + inert disk , (3) 
X + sink -> + sink . (4) 

In practice, we consider a square unit cell of size L x L. In this cell, rid disks 
of unit radius are placed randomly in the following way. The disks are first 
placed periodically, on a square or triangular lattice. A random velocity of unit 
magnitude is assigned to each of them, and the system is evolved as a hard 
ball system during a relaxation time T rdax . The configuration is then frozen 
and n s disks are randomly chosen to be sinks. Each disk has the probability c 
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to be a sink, independently of the others, so n s is always close to end- Three 
examples of such unit cells are shown in Fig. 1. 

The motion of the point particle in the system is defined under periodic bound- 
ary conditions. The collisions on the disks are elastic so that the energy of the 
system is conserved. The magnitude of the point particle velocity is therefore 
a constant, fixed here to one. The particle velocity can be specified by a unique 
coordinate, for example the angle ip between the velocity and the x-axis. Its 
position is specified by two coordinates {x, y}. The three phase-space variables 
of the point particle are therefore T = {x,y,ip}. 



2.2 Escape and fractal repeller 



2.2.1 Asymptotic behavior 

To study the asymptotic time behavior in such a system, we define the survival 
probability P(c,t), which is the probability for a point particle to be still in 
the system at time t, for a sink concentration c now taken to be exactly n s /rid- 
For finite systems in fact this survival probability depends on the positions 
of the scatterers as well as on the distribution of the traps, but as system 
size increases the fluctuations of P(c, t) around its average value will become 
increasingly smaller. Suppose the long time behavior of P(c,t) can be written 
in the form 

P(c,t) ~exp \-Ct a{t) \ . (5) 



In the case of an exponential decay, a would be equal to 1, and for the stretched 
exponential behavior we expect to hold asymptotically in two dimensions, it 
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would be equal to 0.5. This may be generalized into a time dependent a(t) 
defined as [8] 



a(t) 



d 



ln[- In P(c,t)]. 



(6) 



dint 



This quantity, obtained numerically, is plotted versus time in Fig. 2. 

The first graph, in Fig. 2a, represents the values of a for the disk configuration 
represented in Fig. la, but with only 4 traps, corresponding to a trap density 
c ps 0.05. An average has been taken over 10 different trap configurations. 
Figure 2b (respectively 2c) presents the results obtained for the disk configu- 
ration of Fig. lb (respectively lc), with 5 (respectively 6) traps, corresponding 
to c ps 0.05. In all three cases, we observe that, after a transient, a slowly tends 



The escape rate characterizing the exponential decay may be obtained as 
follows [4,5,6]. Let us take iVo initial conditions randomly distributed in the 
phase space of variables T = (x,y,ip) according to some initial measure v . 
After a time t, only a number N t of particles will remain in the system. This 
number will decrease monotonically with t. The set of particles remaining in 
the random Lorentz gas at time t consists of those colliding on a trap, and 
therefore escaping, at an escape time T(r ), which is larger than t: 



The decay of the number of particles is expected to be asymptotically expo- 
nential due to the finite size of the system. The escape rate is thus defined 



to 1 [7]. 



T(f) = {T :t< T(r )}. 



(7) 
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as 

7= lim — In — -= lim — In i/n [T(i)l. (8) 

' t^+oo t N t^+oo t 1 WJ W 

The escape rate is a characteristic quantity of the system, independent of the 
choice of the initial measure u as long as this latter is smooth enough. The 
escape rate is easily accessible by numerical computations, as shown in Fig. 3. 
The logarithm of the survival probability is plotted versus time in Fig. 3 for 
the disk and trap configurations presented in Figs, la and c. The behavior is 
clearly exponential after a short transient. The slopes of the curves give escape 
rates of 7 = 0.025 and 7 = 0.017 respectively. 

2.2.2 Fractal repeller and nonequilibrium probability measure 

A two-dimensional random Lorentz gas of finite spatial extension with traps 
is an open system, obeying an asymptotically exponential decay. Moreover, 
its dynamics is chaotic. It can therefore be studied by the methods devel- 
oped in Refs. [4,5,6], in order to relate its microscopic chaotic dynamics to its 
macroscopic exponential decay. 

An important property of such a system is the existence of a fractal repeller. 
For a typical initial condition r , the particle will collide on a trap and es- 
cape after a finite escape time T(T ). The same will happen if we follow the 
trajectory in the negative time direction. Although a forward and a back- 
ward escape time exist for almost all trajectories, there do exist trajectories 
that never collide on a sink and remain trapped forever in the Lorentz gas in 
both time directions. They can be both periodic and non-periodic. Because 
of the defocusing character of the collisions on the disks, these non-colliding 
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trajectories are unstable and form a fractal set of zero Lebesgue measure in 
phase space: this set is called the fractal repeller T . It is a typical feature of 
chaotic-scattering processes. 

Evidence of the fractal character of this repeller is given by the forward escape 
time T(r ) as a function of the initial condition, as shown in Fig. 4. This 
function is finite for almost all initial conditions since the particle typically 
collides with a sink after a finite time. However, it is infinite for trajectories 
trapped on the fractal repeller. These trajectories start from initial points on 
the stable manifold (in the forward direction) of a trajectory that is trapped in 
both time directions. The infinities of the escape-time function thus establish 
a fractal set formed by the stable manifold of the fractal repeller, W s (J-' n ). 

Fig. 4b-d shows this escape time as a function of the initial position, given by 
an angle 9 around the disk marked by a star in Fig. 4a. The velocity is chosen 
normal to this disk. Due to the random distribution of disks, the fractal set 
is not quasi self-similar, as it was the case in Ref. [6]. However, Figs. 4c-d 
show that a fine structure of escape windows exists at smaller and smaller 
scale. The characteristics of this structure, such as correlation functions and 
the distribution of sizes of empty windows, may be expected to show self- 
similarity on small length scales. 

The fractal repeller is the support of an invariant probability measure which 
can be constructed by considering statistical averages over all the trajectories 
which have not yet escaped during a time interval [— T/2, +T/2] and, there- 
after, taking the limit T — > oo. In this limit the trajectories which do not 
escape during [0, +T/2] have initial conditions on the stable manifold of the 
repeller, W S (J-). On the other hand, the initial conditions of the trajectories 
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which, in the reverse time direction do not escape during [0, — T/2] approach 
the unstable manifold, W^JF). In the limit T — > oo, imposing no escape over 
the whole interval [— T/2, +T/2] selects trajectories which approach closer and 
closer the repeller, given by the intersection W S {!F) fl W^J 7 ) = T . Statistical 
averages over these selected trajectories define an invariant measure having 
the fractal repeller for support. 

In systems with time-reversal symmetric collision dynamics, such as the present 
one, statistical averages calculated over the aforementioned invariant measure 
are equivalent to statistical averages over a conditionally invariant measure 
defined by selecting the trajectories which do not escape during [0,+T/2] 
only and taking the limit T — > oo. This conditionally invariant measure has 
been constructed in Refs. [4,5,6] and we report the reader to these papers for 
further detail. Let us here add that, in the limit of low densities of sinks, the 
system tends to the closed finite Lorentz gas, whereupon the invariant and 
the conditionally invariant measures reduce to the microcanonical equilibrium 
measure. In this limit, the escape process stops and the fractal repeller fills 
the whole phase space. 

2.2.3 Average Lyapunov exponent of the repeller 

The positive Lyapunov of the repeller can be calculated by averaging over the 
N T particles still in the system after time T according to 

X= lim lim - — £ InA T (rJ>), (9) 

where A T at time T is the stretching factor of an initial perturbation on the 
trajectory Tq ^ of surviving particles [4,5,6]. The average (9) is equivalent to an 
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average based on the aforementioned conditional invariant measure as shown 
in Refs. [4,5]. 

Numerically, the positive Lyapunov exponent is obtained as the average slope 
of the logarithm of the stretching factor as a function of time according to Eq. 
(9). This is illustrated in Figs. 5a and 6a, for the disk and trap configurations 
depicted respectively in Figs, la and lc. The average Lyapunov exponent 
respectively equals A = 0.75 and A = 1.67. 

2.2.4 Fractal dimensions of the reseller 

Important characteristics of the fractal properties of a repeller are its gen- 
eralized fractal dimensions. Repellers most often are multifractals, with non- 
trivial generalized fractal dimensions D q [10]. Here the embedding dimension 
of the repeller is three, since we are working in a 3-dimensional phase space. 
Therefore, the fractal dimensions of the repeller must satisfy the inequality 
< D q < 3. Instead of calculating directly the dimensions D q of the fractal 
repeller of the full flow we will consider a set of so-called partial fractal dimen- 
sions d q defined for the intersection of the stable manifold of the repeller T 
with a line C across this stable manifold. Like T this mostly is a multifractal 
set. Obviously the values of d q must be between and 1. Because the system is 
time-reversal symmetric the partial dimension in the stable direction is equal 
to the one in the unstable direction. The direction of the flow contributes a 
partial dimension equal to unity. Hence for given q, the dimension D q is related 
to the partial dimension d q by D q = 2d q + 1. 

In the present paper, we focus on the fractal dimensions for q = and q — 1. 
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For q = 0, we have the Hausdorff dimension commonly defined by 

d H = d = lim — ^ , 10 

£^o In e 

where N{e) is the number of intervals Ij of equal size e, required to cover 
the fractal / defined on a line. For q — 1, we find the so-called information 
dimension 

I Ni 

di = d x = lim — Y,PA £ ) ln Pj( e ) , (11) 

where Pj(e) = /i(Ij fl /) is the probability weight of the interval Ij of size e. 
Here, /i denotes the probability measure describing the distribution of points 
in the fractal set /. The Hausdorff and information codimensions are defined 
respectively as c# = 1 — d H and cj — \ — dj. 

For fractal repellers of an open hyperbolic system, a relation exists between the 
partial information codimension cj, the (largest) positive Lyapunov exponent 
A, and the escape rate 7 [4,5,6] 

7 = A cj. (12) 



To a good approximation the information codimension can be replaced by the 
Hausdorff codimension [4,5]. This last dimension may be calculated with the 
aid of a numerical algorithm developed by the group of Maryland [4,5,11]. 
The basic idea of this algorithm is to consider an ensemble of pairs of initial 
conditions, separated by e, along the line £ defined above. A pair is said to 
be uncertain if there is a discontinuity of the escape-time function between 
both initial conditions. On the other hand, when the pair is certain, the initial 
conditions belong to an interval of continuity of this function. An estimation of 
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the fraction of uncertain pairs may be obtained in the following way: typically 
for an uncertain pair the time after which the pair becomes uncertain coincides 
with the time that the distance between the points becomes of order unity. 
This time may thus be estimated as t u = (1/A) ln(l/e). Since trajectories 
escape at the rate 7, the fraction of uncertain pairs would be proportional to 
exp(— jt u ) ~ e 7 / A = e Cl . However, such an estimation needs to be refined to 
take into account the multifractal character of the repeller and the possibility 
that Cf 7^ ch due to local fluctuations of the Lyapunov exponents. A refined 
estimation of the fraction of uncertain pairs given in the Appendix A for the 
simple case of a piecewise-linear one-dimensional map shows that 

f(e)~e c «, (13) 



as proved for the general case in Ref. [11]. We notice that, to leading or- 
der in the sink density c, the Hausdorff codimension equals the information 
codimension [6]. 

In the case of the Lorentz gas, discontinuities in the escape-time function 
are not very practical to calculate numerically. Instead one may calculate the 
fraction of pairs that do not collide with exactly the same sequence of disks 
until escape from the system. For small e this will behave in precisely the same 
way as the fraction of uncertain pairs, because it is determined by the same 
criterion of reaching a distance of order unity between both trajectories of a 
pair before escaping. The details of the algorithm used here are described in 
Refs. [4,5]. 

In Figs. 5b and 6b, the logarithm of the fraction of uncertain pairs (in fact 
pairs with different collision sequences before escape) is plotted as a function of 
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the logarithm of the small e separating the two initial conditions of a pair, for 
the systems defined respectively by the unit cells depicted in Figs la and lc. 
After a transient regime, the linear behavior confirms the power-law prediction 
of Eq. (13) . The slope gives a value of the Hausdorff codimension equal to 
ch = 0.034 and ch = 0.01 respectively. 

For both configurations, the value of the product \c H is in good agreement 
with the value of the corresponding escape rate: Xch = 0.025 = 7 in the 
former case and Xch = 0.017 = 7 in the latter. Since ch = cj to leading order 
in the sink density c, this confirms the validity of the relation (12) for random 
Lorentz gases of finite spatial extension. This expression connects directly two 
quantities characterizing the microscopic chaotic dynamics of the Lorentz gas, 
the average Lyapunov exponent A and the Hausdorff codimension ch of the 
fractal repeller, and a quantity describing the macroscopic transport behavior 
of this open system, the escape rate 7. 

3 The two-dimensional random Lorentz gas of infinite spatial ex- 
tension 

3. 1 Description of the model and simulation method 

We now consider a two-dimensional random Lorentz gas of infinite spatial 
extension. This supposes that the sinks are randomly distributed among an 
infinite number of disks. As a consequence, the density of larger and larger 
trap-free regions decays exponentially with their volume but does not vanish. 
A suitable method is needed to simulate this situation. It cannot be done with 
the simple method presented in the previous section: if, in a possibly large, 
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but in practice always finite system, the nature of the disks is fixed a priori, 
the size of the largest trap-free region is determined. Trap-free regions of a size 
that is very rare for the given system size in practice just will not be observed. 
However, a certain probability c to be a sink can be assigned to each disk [8,7]. 
This leads to an average over the different possible sink configurations, as a 
result of which the probabilities for finding large trap-free regions will attain 
just the proper values. 

Let us define the corresponding configuration-averaged survival probability 
P(c, t), where c may be interpreted alternatively as the average concentration 
of sinks. The time evolution of P(c, t) is given by 

P(c, t) = E(l - c) n P(n, t) = ((l- cD (14) 

n 

where n = n t is the number of different disks on which the point particle has 
collided during the time t and p(n, t) is the probability that this number is 
equal to n. (...) is an ensemble average. (1 — c) is the probability that a disk 
is not a sink, so that (1 — c) n is the probability that none of the disks met at 
time t is a sink, in other words the probability that the point particle survives. 

To compute the probability p(n, t) to meet n different disks in a time t, a 
purely deterministic simulation is used. The two-dimensional random Lorentz 
gas in which the point particle diffuses is defined as in the previous section, 
by a square unit cell of size L x L. Periodic boundary conditions are applied. 
However, an index is associated with each square cell, allowing us to know 
which one is currently visited. This is necessary to distinguish one disk from 
its copies in the other cells. Indeed, if the point particle meets a disk and, 
later, one of its copies, this should be considered as collisions with different 
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disks. In other words, although the scatterers are located on a periodic lattice, 
be it one with a very large unit cell of random filling, the distribution of traps 
has no periodicity whatsoever. 

3.2 Theoretical predictions and simulation results 

The asymptotic time dependence of the survival probability P(c, t) is expected 
to be of the form (1) [2,3,8,12]. For checking this dependence in our system, 
the easiest method using the idea of Eq. (14) is a direct simulation, in which 
p(n, t) is estimated as an average over a large number of uniformly distributed 
initial conditions. The system considered here has as unit cell the one depicted 
in Fig. lc, of size 25 x 25 and with 127 disks. The number of initial conditions 
is of the order of 10 s . 

We present our data following a procedure proposed by Barkema et al. in Ref. 
[12]. At short times and small concentrations, Eq. (14) can be approximated 
by [13] 

P(c,t) = ((l-c) rit ) ~ (1 -c) {nt) . (15) 

This is the so-called Rosenstock approximation. In two dimensions, (n) scales 
as AnxpDtj hit, with p the average density of scatterers. We therefore expect 
a time dependence of the survival probability of the form: 

— In [P(c, t)] ~ AirpDut/ lnt at short times and small concentrations, 

- In [P(c, t)] ~ (ut) 1/2 asymptotically, (16) 

with again u — — ln(l — c). An efficient way to visualize the crossover from 
one regime to the other is to plot — In [P(c, t)] / lnt as a function of y/ui/ lnt 
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[12,14]. In the Rosenstock regime, we expect a quadratic dependence, while 
the stretched exponential regime should produce a linear behavior. The results 
are represented in Fig. 7 in a double- logarithmic plot. The data collapse for 
values of c taken in the range [0.01, . . . , 0.09, 0.1, . . . 0.9]. The crossover from 
a slope 2 to a slope 1 is clear. 



The intersection point of the two linear fits gives an estimation of the crossover 
time from the Rosenstock regime to another regime which ends with the 
Donsker-Varadhan asymptotics. Using this as reference point, we obtain the 
dependence on the trap density c of the crossover time and of the crossover 
survival probability [12], as presented in Figs. 8a and b respectively. These re- 
sults show that the crossover occurs at shorter and shorter times and at larger 
and larger probabilities as the trap density increases. We observe that the 
crossover probability changes from O(10~~ 20 ) to (9(10~ 5 ) as the trap density 
only varies from 0(1O- 2 ) to 0(1). 

A comment is here in order. The extensive numerical simulations of Ref. [15] 
show that the Donsker-Varadhan behavior is approached only very gradually. 
It is known to provide an exact asymptote for the logarithm of P(c, t) including 
the prefactor of (ut) 1 ^ 2 . If one inserts the relevant values of the parameters 
one find this asymptote in Fig. 7 about 20% above the asymptote obtained 
from the data fit. Given the limited time range over which our numerical data 
were obtained we may conclude that their behavior is fully compatible with 
the one observed for stochastic systems in Ref. [15]. 
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4 Transition from finite to infinite random systems 

We have considered finite random Lorentz gases with exponential behavior 
and infinite ones with stretched exponential behavior. How does the transition 
between these two regimes evolve with increasing system size LI And what 
happens to the fractal repeller? 

To answer these questions, let us first briefly sketch the arguments of Grass- 
berger and Procaccia for the infinite case [3]. They consider diffusing particles 
distributed according to a density p c (r,t). The average density n c (r,t) is de- 
fined as 

n c (r,t) = (p c (T,t)) a (17) 

where (...) s is an average over all the realizations of trap distributions or 
equivalently over all choices of the origin from which all particles start diffusing 
initially. A perfectly absorbing sphere of large radius R centered at the origin 
is added to the system. The value of the new density n' c (r,t) obtained in the 
presence of this sphere is necessarily smaller than n c (r,t) 

n c (r,t)>n' c (r,t). (18) 

n' c {r,t) depends on the distribution of traps inside the sphere. If {S} = 
{Si, S 2 , ■ ■ ■ , Sk} are the positions of the traps in the sphere, we have 

n c (v,t) > n' c (v,t) = XX(r,f|{S}) P({S}). (19) 

{S} 

The sum of the right-hand side can be bounded from below by considering 
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only one term, corresponding to the situation of a sphere free of traps: 

n c (r,t)>ri e (r,t\{S} = 0)P(0). (20) 



A random distribution of sinks corresponds to a Poisson distribution. The 
probability of finding k traps in a volume v is given by 

N k 

P k = —exp(-N) (21) 



with average value N = n s v, where n s is the average density of sinks. The prob- 
ability to have no trap in the sphere is thus given by P(0) = exp(— C n s R d ), 
where d is the dimension of the system and C the volume of the unit ball in 
d dimensions. 

On the other hand, n' c (r,t\{S} = 0) is obtained as solution of a diffusion 
equation of coefficient V with absorbing boundary conditions on the border 
of a sphere of radius |r| = R. Such a solution is known to decay asymptotically 
as the exponential exp(— Kdj^t), where Kd is a constant depending on the 
dimensionality of the system. 

Therefore, Eq. (20) becomes 

n c (r,t) > n' c (r,t\{S} = 0)P(0) ~ exp (- Kd ^ t - Cn s R d ^ . (22) 



The right-hand side may be maximized as a function of R, with R between 
and +oo in an infinite system. This gives us an optimal value for R 

R= 2Kd^> t \ =at l/(*-2)_ (23) 

V dC n s J 
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Inserting this result in Eq. (22), one obtains 

n c (r, f) > exp [-[3t d/{d+2) ] . (24) 

Let us now consider a finite system of linear size 2L with a fixed configuration 
of traps. As long as R = at l/d+2 < (^) , the optimal value of R is found 
in a fraction of the volume described by (21) and the stretched exponential 
behavior is observed. However, for longer times, the maximal value of R is 
proportional to l rs (In Ljn s ) l l d . Inserting R — l in Eq. (22) leads to an 
exponential decay of the form exp (^—K d ^ tj. As L increases, the time t* at 
which the transition occurs increases. 

A logarithmic dependence of t* on L in two dimensions has been proven in 
Ref. [7]. 

We now want to study the properties of the repeller. As a preparation, let us 
first consider the simple one-dimensional map with escape represented in Fig. 
9. Starting from the full unit interval at time zero one finds after k discrete 
time steps an approximate forward repeller consisting of Nk = 2 k pieces of 
sizes Ek = (l/3) fc . The number of pieces is equal to the exponential of the 
Kolmogorov-Sinai entropy times time, N k = exp(tiKsk), and the linear size of 
these pieces is given by Ek = exp(— Xk) with the Lyapunov exponent A = In 3. 
This last equality can be intuitively understood as follows: as time goes on, new 
structures appear on the repeller on smaller and smaller scales. Therefore, long 
times correspond to small spatial scales. The total length Lk of the repeller at 
time step k is the product of Nk and Ek- It decreases exponentially with time 
as 

L k = N k E k = exp(h K sk) exp(-Afc) = exp(-7fc), (25) 
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where 7 = A — hxs is the escape rate. 



The fractal (Hausdorff) dimension of the repeller is thus given by the well- 
known result 



Dh = lini 



. h\N(e) h KS In 2 



oin(i) 



A ln3' 



(26) 



in the illustrative example of Fig. 9. On the other hand, we know that the 
fraction of uncertain pairs f(e), as defined for the Maryland algorithm, is 
proportional to the size e to the power C#, with Ch the codimension of the 
repeller. From this we obtain 



f{e) ~ e Ch ~ e l -° H ~ ~ eN{e). 



(27) 



Let us now imagine a one-dimensional system in which the escape would obey 
a stretched exponential law. In this case, the remaining length of the repeller 
at time k would obey 



L k = N k e k = exp(-c k 1/2 ). 



(28) 



iFrom the relation e = exp(—k A), we can express the time as k = ^ In K We 



get 



L k exp(-c k 1/2 ) 1 
N(e) = — = = - exp 

see 



• [\ ln-^ 

A e 



(29) 



/,From Eq. (27), we have 
f(e) ~ eN(e) ~ exp 



(30) 



This dependence of the fraction of uncertain pairs on e for a system with 
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stretched exponential behavior can be tested numerically for the two-dimensio- 
nal infinite random Lorentz gas, with the help of the Maryland algorithm. 
The algorithm has to be slighty modified to simulate an infinite system. As 
mentioned before, the absorbing or non-absorbing nature of the disks is not 
fixed a priori, but each disk has a probability c to be a sink. At each collision on 
a disk, the separation of the pair of particles considered is tested. If separation 
occurs, the fraction of uncertain pairs is increased by (1 — c) n , where n is the 
number of different disks that the pair has met before separating, so (1 — 
c) n is the probability that the particles have not escaped before separating. 
Simultaneously, the fraction of escaped pairs is increased by [1 — (1 — c) n \. 
In Fig. 10a, the logarithm of the fraction of uncertain pairs, log 10 /(e), is 
plotted as a function of the logarithm of the distance between the two initial 
conditions of the pair, log 10 e. In order to check that this curve can be fitted by 
the expression (30), we plot in Fig. 10b log[— log 10 /(e)] versus log(— log 10 e). 
We obtain a linear behavior with slope |, in excellent agreement with Eq. 
(30). 

In the case of a system of finite spatial extension but with high density of 
disks, as represented in Fig. lc, we see in Fig. 3b that the transient behavior 
persists until t ~ 50. According to the relation e = exp(tA), this corresponds 
to log 10 e ~ —34. This is consistent with the deviation from the linear curve 
which is observed in Fig. 6b for large values of e. 

In two-dimensional systems like the Lorentz gases we consider here, the codi- 
mension appearing in the relation /(e) ~ e CH is a partial codimension, asso- 
ciated with a partial dimension dn- The total dimension is given by D H = 
1 + 2dn- The total number of small three-dimensional pieces needed to cover 
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the fractal repeller is given by 



N(e) ~ - 

£ 



2 i r n i\ 1/2 

~ — exp 



e 3 



-2c ( - In 

,A e 



(31) 



This type of dependence of N(e) on the size e leads to an integer value of the 
Hausdorff dimension D H . Indeed, 



2 



D h = ^^t^ = ^-^c\-^- t _) =3. (32) 



+ ° ln ft) 



This dimension does not contain the information concerning the correction to 
the power law dependence appearing in Eq. (31). We notice that corrections 
to a power law also appear for Brownian motion in the plane [16] . In this case, 
the dependence essentially would take the form 

N ^ ~ peei < 33) 

e 

leading to the integer Hausdorff dimension D H = 2. 

When the size L of a finite system with fixed traps increases, the time during 
which stretched exponential behavior is observed increases likewise, and the 
range of e over which the fraction of uncertain pairs deviates from a power-law 
behavior in e. Correspondingly, as seen from Fig. 10, the slope in the linear 
regime of log 10 f(e) versus log 10 e becomes smaller and smaller, i.e., the partial 
Hausdorff codimension of the repeller will be smaller. This means that with 
increasing L the Hausdorff dimension of the repeller progressively increases 
towards the asymptotic value Dh = 3. The same holds for the information 
dimension, as in the dominant regions the trap density approaches zero and 
we saw before that both dimensions approach each other in this limit. 

In the case of the two-dimensional random Lorentz gas of infinite spatial ex- 
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tension, the integer value Dh = 3 of the dimension of the fractal repeller can 
be understood by intuitive arguments. Indeed, this dimension corresponds to 
the case of a closed system, without escape. In an infinite system, for a given 
time t, there exists with a finite probability a region free of traps of size 
R = [-^Q^f-tj . The point particles in this region behave as particles in 
a closed system until they escape through the "boundary" at R. As time goes 
on, larger and larger trap-free regions have to be considered. Even if their 
density decays exponentially with their size, it does not vanish, so that these 
regions give the dominant contribution to the measure of the repeller. 

Our results show that the repeller is indeed a fractal object in spite of its 
integer dimension D H = 3. This fractality is characterized by the complicated 
scaling behavior of the fraction (30) of uncertain pairs or of the number (31) 
of pieces needed to cover the fractal. These theoretical predictions are nicely 
confirmed by our numerical results. 

5 Conclusions 

We have studied in this paper several types of two-dimensional random Lorentz 
gases with sinks. 

The first model considered is a two-dimensional random Lorentz gas with 
sinks, of finite spatial extension. In this system, the sinks are randomly dis- 
tributed among the disks. Due to the finite size of the system, the size of 
the largest trap-free region is bounded. Therefore, the decay of the number of 
particles is asymptotically exponential. For this open system, we have shown 
the existence of a fractal repeller formed by the trajectories that never escape 
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from the system. Using a non-equilibrium measure on the repeller, we have 
characterized its chaotic and fractal properties. We have obtained the escape 
rate 7, the positive Lyapunov exponent A, and the fractal dimensions of the 
repeller. A relation between 7, A, and the Hausdorff codimension ch of the 
fractal repeller has been numerically tested in the limit of low sink densi- 
ties. This relation establishes a link between the macroscopic decay and the 
underlying microscopic chaos. 

The second case considered is the two-dimensional random Lorentz gas with 
sinks, of infinite spatial extension. This situation can be simulated by giving 
to each disk a certain probability to be a sink, while distinguishing between 
different periodic images of the same disk. This leads to an average over the 
different sink configurations. The decay of the survival probability asymptot- 
ically obeys the stretched exponential law (1). Here this behavior has been 
studied numerically for the first time in a deterministic chaotic system. 

For a dense 2D random Lorentz gas of finite size L, the transient stretched 
exponential regime persists during a time fixed in first approximation by the 
typical size of the largest trap-free region in the system [see Eq. (23)]. As L 
increases, the transient regime is observed for longer and longer times. Simul- 
taneously, the dimension of the repeller increases, reaching the integer value 
D H = 3 in the infinite-size limit. We have here been able to show that the 
repeller remains fractal in the limit and to answer the puzzle mentioned in the 
introduction by finding the nontrivial scaling property of this plain fractal. Its 
fractal property can be characterized equivalently by the fraction of uncertain 
pairs of trajectories starting from initial conditions separated by e or by the 
number of cells of size e needed to cover the fractal. Both scaling functions 
show a correction to a simple power law proving that the repeller remains 
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fractal. The correction depends on the Lyapunov exponent as well as on the 
exponent 1/2 of the stretched exponential law (1) in d = 2. This nontrivial 
scaling is remarkably confirmed by our numerical results. 

In conclusion, we have been able to show that, in a random deterministic 
chaotic system of infinite spatial extension with sinks, the repeller composed 
of all the trajectories that never meet a sink is a plain fractal with special 
scaling properties related to the stretched exponential decay of the survival 
probability. 
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A Direct proof that the uncertainty exponent is the Hausdorff 
codimension 

Let us consider the piecewise-linear one-dimensional map: 



[Ai(l-x n ) A<x n <l 

We assign the symbol oo = to the interval [0, A] and the symbol uo = 1 to 
the interval [A, 1]. 




< x n < A 



(A.l) 
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If A , Ai > 2, almost all the trajectories escape after a finite number of itera- 
tions. The subinterval I CJ0U ; 1 .. Mn _ 1 is defined as the set of initial conditions for 
the trajectories which escape in n iterations and successively visit the intervals 
corresponding to the symbolic sequence ujQUJi...uj n -\. These subintervals form 
a partition of the unit interval [0, 1]: 

Ua; u;i...u; n -i Iu> Q u>i...u> n _ 1 [0, 1] (A-2) 

The length of the subinterval 1^1...^-! i s 

K 

^woUi...u„_i = "7 7 7 (A-3) 

Since the sum of lengths should be equal to one 

oo 

E E ^WoWi—Wn-i = 1 (^-4) 

n=l W0 u l- W n-1 



the constant fT is given by 

K = ^y/- 1 (A.5) 
A 1 + Ai 1 



and the previously introduced constant A can now be fixed to the value: 

A = I+* (A.6) 
Jin 



A pair is uncertain if it falls in a subinterval which is smaller than e. Therefore, 
the fraction of uncertain pairs is equal to 

oo 

/(e) = E E W.^-l ■ (A.7) 

„=lu;oa;i...a;„_i ^o<n-"n-i <£ 

The sum can be reordered in terms of the numbers no (resp. n\) of symbol 
(resp. 1) in the symbolic sequence ujQUJi...uj n - 1 . We have that n + ri\ = n. 
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The fraction of uncertain pairs becomes 
^ y, (np + m)! K 

„ =0n 1= n l n l ! A A l 



f(r) = V V K + K ( A 



We use the steepest-descent method for this double sum. The domain of sum- 
mation is the set of pairs of non-negative integers (n ,ni) above the straight 
line nolnA + nilnAi = ha(K/e) in the plane (n ,ni). The function to be 
summed decreases exponentially to zero as the exponents n and n\ tend to 
infinity above the straight line. Therefore, the sum will be dominated by the 
maximum of the function on the straight line. If we introduce the alternative 
variables 



we thus find 



/(e)~e(-l -e l - D (A.10) 



with 



D = x m h^iln (A. 11) 

x Xi 



and the constraint 



z lnAo + ZilnAi = 1. ( A - 12 ) 



The maximum of the function is given by the steepest-descent condition 

dD . x + xi dx 1 x + x 1 ,.„, 
- — = In h - — In = 0. (A. 13) 

ax Q x ax xi 
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With the constraint (A. 12), we find that the steepest-descent condition reduces 
to 

In Ax ln^±^ = lnA ln^°±^. (A.14) 

x x 1 

We now show that the exponent D is equal to the Hausdorff dimension. In 
the present piecewise-linear map, the Hausdorff dimension of the repeller must 
satisfy by definition the following condition 

aTxt 1 - (A ' 15) 

Let us show that it is indeed the case. For this purpose, we first express the 
quantity D in terms of x and x\ by Eq. (A. 11). Then, we use the steepest- 
descent condition (A.14) and, thereafter, the constraint (A. 12). The calcula- 
tions are the following: 



J_ J_ _ -x in ^±£1 lnAo _ xl ln ad^l lnAo 

A? Af 

+e *i (A. 16) 

= e -(x lnA +^ilnA 1 )ln^fl±£i + ^(xolnAo+xi InAi) ln ^ ^ 

_lnfH±£l _in£n±£l , , 

= e + e =i=l (A. 18) 



Q.E.D. 

The previous calculation provides a direct proof that the uncertainty exponent 
of the Maryland algorithm is the Hausdorff codimension in a piecewise-linear 
map with two branches. Consequently, the uncertainty exponent should be 
expected to differ from the information codimension on general multifractal 
repellers. Indeed, by using the example (A.l), we have numerically observed 
that the uncertainty exponent agrees with the Hausdorff codimension but 
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differs from the information codimension if Ao ^ Ai. 

The direct proof presented here can be extended to piecewise-linear one- 
dimensional maps with more than two branches. In the case of to branches, 
there would be to variables rii or Xi and (to — 1) steepest-descent conditions 
such as (A. 13) because of the constraint that the maximum point belongs to 
the hyperplane Y^Li fiilnAj = \n(K/e). Otherwise, a Lagrange multiplier may 
be introduced, to deal with this constraint in the search of the maximum. In 
the general case, the argument of Ref. [11] provides the proof. 
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(a) 




(b) 



S3 * ss^Shdf-*' 



(c) 




Fig. 1. Typical unit cell of a finite two-dimensional random Lorentz gas: (a) L = 25, 
n d = 72, n s = 9; (b) L = 25, n d = 88, n s = 7; (c) L = 25, ra d = 127, n s = 23. The 
filled disks are the sinks. 
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(c) 

Fig. 2. Two-dimensional random Lorentz gas of finite spatial extension; Effective 
exponent a plotted versus time, for a trap density c as 0.05: (a) in the case of the 
disk configuration represented in Fig. la, but with n s = 4 traps; (b) in the case of 
the disk configuration represented in Fig. lb, but with n s = 5 traps; (c) in the case 
of the disk configuration represented in Fig. lc, but with n s = 6 traps. 
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(a) (b) 

Fig. 3. Two-dimensional random Lorentz gas of finite spatial extension; Logarithm 
of the survival probability as a function of time: (a) for the disk and trap config- 
uration represented in Fig. la, (b) for the disk and trap configuration represented 
in Fig. lc. After a transient time the behavior is exponential and the slope of the 
curve gives the value of the escape rate 7. 
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(d) 



Fig. 4. Escape time as a function of the initial position: (a) The initial positions are 
taken at a varying angle 8 around the disk marked by a star; the initial velocity 
is always normal to the disk. Escape time for (b) 9 G [0, 2ir]; (c) G [27r/3,7r]; (d) 
9 G [57r/6,7r]. 
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(a) (b) 

Fig. 5. Two-dimensional random Lorentz gas of finite spatial extension denned by 
the unit cell represented in Fig. la: (a) Average of the logarithm of the stretching 
factors, as a function of time t. The slope gives the average Lyapunov exponent 
A = 0.75. (b) Maryland algorithm: the logarithm of the fraction of uncertain pairs 
as a function of the logarithm of the difference e between the two initial conditions of 
the same pair. The slope gives the Hausdorff codimension ch = 0.034. The product 
Xch = 0.025 is in good agreement with the value of 7, obtained from Fig. 3a. 
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Fig. 6. Two-dimensional random Lorentz gas of finite spatial extension denned by 
the unit cell represented in Fig. lc: (a) Average of the logarithm of the stretching 
factors, as a function of time t. The slope gives the average Lyapunov exponent 
A = 1.67. (b) Maryland algorithm: the logarithm of the fraction of uncertain pairs 
as a function of the logarithm of the difference e between the two initial conditions 
of the same pair. The slope gives the Hausdorff codimension ch = 0.01. The product 
\ch = 0.017 is in good agreement with the value of 7, obtained from Fig. 3b. 
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(ut) /In t 

Fig. 7. Infinite two-dimensional random Lorentz gas: Crossover from the Rosenstock 
regime to the stretched exponential regime. — In [P(c, t)} / In t is plotted as a function 
of \fwtj In t in a double logarithmic plot. The solid lines are fits to the data with 
slopes 2 and 1. They cross at the point (1.7,2.7). 




(a) (b) 

Fig. 8. (a) Crossover time from the Rosenstock regime to the stretched exponential 
regime as a function of the trap density c. (b) Dependence of the crossover survival 
probability on the trap density c. 
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(b) 

Fig. 10. Maryland algorithm in the case of a random Lorentz gas of infinite spa- 
tial extension: the system considered is represented in Fig. lb, with a sink density 
c = 0.6. (a) The logarithm of the fraction of uncertain pairs is plotted as a func- 
tion of the logarithm of the difference e between the two initial conditions of the 
same pair, (b) Logarithmic representation of the same quantities, confirming the 
dependence predicted by Eq. (30). 



40 



